Data Analysis Notebook

#Authors: Ian McElveen and Cecilia Gonzales
#Author Date: 6/9/25

#Purpose: The purpose of this notebook is to house all data set transformation, cleansing, visualization, statistical analysis, and note-taking for the 2025 CU Athletic Department Sports Science Internship Program

#LAST UPDATED: 6/12/2025

#Including helpful libraries
library(tidyverse)
library(readxl)
library(aod)
library(dplyr)
library(ggplot2)
library(lubridate)
library(boot)
library(ROCR)
library(purrr)
library(ggforce)
library(lme4)
library(cv)
library(caret)
library(rsample)
library(yardstick)
library(knitr)


#Example of importing our first data set from a CSV file to a data frame.

# df <- read.csv("./foo/bar.csv")

#Load in MBB - Statistic Tracking Report csv
stats <- read.csv("data-sets/data-sets-uncompressed/data-sets-compressed/Reactive Strength Index/MBB - Statistic Tracking Report.csv")

# Load in ACWR - Kinexon MBB
kinexon_ACWR <- read.csv("data-sets/data-sets-uncompressed/data-sets-compressed/Reactive Strength Index/ACWR - Kinexon MBB.csv")

#Loading in VALD - Dynamo MBB
VALD_dynamo <- read.csv("data-sets/data-sets-uncompressed/data-sets-compressed/Reactive Strength Index/VALD - Dynamo MBB.csv")

# Load in VALD - ForceDecks MBB csv
VALD_ForceDecks <- read.csv("data-sets/data-sets-uncompressed/data-sets-compressed/Reactive Strength Index/VALD - ForceDecks MBB.csv")

# Load in Kinexon Session
kinexon_session <- read.csv("data-sets/data-sets-uncompressed/data-sets-compressed/Reactive Strength Index/Kinexon Session MBB.csv")

Data Cleaning

# Remove columns where all data points are NA
stats_clean <- stats |>
  select(where(~ !all(is.na(.))))


# If all three date columns are equal, remove the redundant ones
if (all(stats_clean$Date == stats_clean$Date.Reverse, na.rm = TRUE) &&
    all(stats_clean$Date == stats_clean$Session.Date, na.rm = TRUE)) {
  
  stats_clean <- stats_clean |>
    select(-c(Date.Reverse, Session.Date))
}


# Convert if Date column is NOT already in Date format
if (!inherits(stats_clean$Date, "Date")) {
  stats_clean <- stats_clean %>%
    mutate(Date = as.Date(Date, format = "%m/%d/%Y"))
}
# Confirm it's now in Date format
class(stats_clean$Date)

[1] "Date"

stats_clean <- stats_clean[!is.na(stats_clean$Plus..Minus),]

# Add more stats
stats_clean <- stats_clean %>%
  filter(Date >= as.Date("2024-09-20", '%Y-%m-%d')) %>%
  group_by(anon_id) %>%
  mutate(Player.Average.Rebounds = mean(na.omit(Rebounds)),
         Player.Average.Steals = mean(na.omit(Steals)),
         Player.Average.Assists = mean(na.omit(Assists)),
         Player.Average.Turnovers = mean(na.omit(Turnovers)),
         Player.Average.Points.Scored = mean(na.omit(Points.Scored)),
         Player.Median.Rebounds = median(na.omit(Rebounds)),
         Player.Median.Steals = median(na.omit(Steals)),
         Player.Median.Assists = median(na.omit(Assists)),
         Player.Median.Turnovers = median(na.omit(Turnovers)),
         Player.Median.Points.Scored = median(na.omit(Points.Scored))) %>%
  ungroup() %>%
  group_by(Date) %>%
  mutate(Team.Rebounds = sum(na.omit(Rebounds)),
         Team.Steals = sum(na.omit(Steals)),
         Team.Assists = sum(na.omit(Assists)),
         Team.Turnovers = sum(na.omit(Turnovers)),
         Team.Points.Scored = as.numeric(str_split(Game.Score, "-")[[1]][1])) %>%
  ungroup() %>%
  mutate(Team.Average.Rebounds = mean(Team.Rebounds),
         Team.Average.Steals = mean(Team.Steals),
         Team.Average.Assists = mean(Team.Assists),
         Team.Average.Turnovers = mean(Team.Turnovers),
         Team.Average.Points.Scored = mean(Team.Points.Scored)) %>%
  mutate(Player.Good.Game = ifelse(Rebounds >= Player.Median.Rebounds & Steals >= Player.Median.Steals & Assists >= Player.Median.Assists,1,0),
         Player.Bad.Game = ifelse(Rebounds < Player.Median.Rebounds & Steals < Player.Median.Steals & Assists < Player.Median.Assists,1,0),
         Team.Good.Game = ifelse(Team.Rebounds >= median(Team.Rebounds) & Team.Steals >= median(Team.Steals) & Team.Assists >= median(Team.Assists),1,0),
         Team.Bad.Game = ifelse(Team.Rebounds < median(Team.Rebounds) & Team.Steals < median(Team.Steals) & Team.Assists < median(Team.Assists),1,0)) %>%
  select(-c(Sport, Day.of.Week, Month, Month.with.Number, Year, Number.of.Competitions.Won, Total.Daily.Competitions.Available, Total.Competitions.Won, Total.Game.Minutes, Offensive.Rebounds, Defensive.Rebounds, Buff.Plays, Percentage.Daily.of.Competitions.Won)) %>%
  filter(!(anon_id %in% c("ID_10", "ID_14", "ID_15", "ID_18", "ID_54", "ID_64")))

# Remove columns where all data points are NA
kinexon_ACWR_clean <- kinexon_ACWR |>
  select(where(~ !all(is.na(.))))



# If all three date columns are equal, remove the redundant ones
if (all(kinexon_ACWR_clean$Date == kinexon_ACWR_clean$Date.Reverse, na.rm = TRUE) &&
    all(kinexon_ACWR_clean$Date == kinexon_ACWR_clean$Session.Date, na.rm = TRUE)) {
  
  kinexon_ACWR_clean <- kinexon_ACWR_clean |>
    select(-c(Date.Reverse, Session.Date))
}


# Convert if Date column is NOT already in Date format
if (!inherits(kinexon_ACWR_clean$Date, "Date")) {
  kinexon_ACWR_clean <- kinexon_ACWR_clean %>%
    mutate(Date = as.Date(Date, format = "%m/%d/%Y"))
}
# Confirm it's now in Date format
class(kinexon_ACWR_clean$Date)

[1] "Date"

# Selecting only relevant columns
kinexon_ACWR_clean <- kinexon_ACWR_clean %>%
  select(anon_id, Date, Day.of.the.Week, Month, Year, Position, Session.Type, Session.Description, Practice.Classification, Last.7.Day.Acceleration.Load, Today.s.Acceleration.Load, Acceleration.Load.ACWR, Acceleration.Load.EWMA.ACWR, Acute.Total.Acceleration.Load, Acute.Total.Acceleration.Load.EWMA, Chronic.Total.Acceleration.Load, Chronic.Total.Acceleration.Load.EWMA, Total.Low.Acceleration.Load, Total.Medium.Acceleration.Load, Total.High.Acceleration.Load, Total.Very.High.Acceleration.Load, Last.7.Day.Total.Jump.Load, Today.s.Jump.Load, Total.Jump.Load.ACWR, Acute.Total.Jump.Load, Acute.Total.Jump.Load.EWMA, Chronic.Total.Jump.Load, Chronic.Total.Jump.Load.EWMA, Last.7.Day.Total.Minutes, Today.s.Minutes, Minutes.ACWR, Minutes.EWMA.ACWR, Acute.Minutes, Acute.Minutes.EWMA, Chronic.Minutes, Chronic.Minutes.EWMA, Total.Distance.EWMA.ACWR)

#selecting important variables, removing rows with missing data, only selecting hand movements
VALD_dynamo_clean <- VALD_dynamo %>%
  filter(Body.Region == "Hand") %>%
  mutate(Right.Side = ifelse(Repetition.Type.Laterality=="RightSide",1,0)) %>%
  select(anon_id,   Date, Max.Force.Newtons, Avg.Force.Newtons, Max.Impulse.Newton.Seconds,  Max.Rate.Of.Force.Development.Newtons.Per.Second,  Avg.Rate.Of.Force.Development.Newtons.Per.Second,   Min.Time.To.Peak.Force.Seconds, Start.Offset.Seconds, Repetition.Duration.Seconds,  Repetition.Max.Force.Newtons,   Impulse.Newton.Seconds, Rate.Of.Force.Development.Newtons.Per.Second,   Time.To.Peak.Force.Seconds,Right.Side) %>%
  na.omit()

# Convert if Date column is NOT already in Date format
if (!inherits(VALD_dynamo_clean$Date, "Date")) {
  VALD_dynamo_clean<- VALD_dynamo_clean%>%
    mutate(Date = as.Date(Date, format = "%m/%d/%Y"))
}
# Confirm it's now in Date format
class(VALD_dynamo_clean$Date)

[1] "Date"

#This data set contains information from each individual hand in a strength test meaning that there are two measurements per day for each player. This makes it hard to understand in plots especially when plotted like time series data. For future reference, it might be useful to take an average of both sides or note that there is a difference in each hand and plot the left and rights sides separately.

#### Note that you have to group on id and the date in order to get a valid average of both sides for a specific day for each individual player. 

VALD_ForceDecks_clean <- VALD_ForceDecks %>%
  select(where(~ !all(is.na(.)))) %>%
  mutate(X=1:3384) %>%
  filter(Test.Type == "DJ") %>%
  select(X,anon_id, Date, Week, Position, Test.Type, Weight..kg., Trial, Flight.Time..s., Jump.Height..Flight.Time...cm., Eccentric.Concentric.Mean.Force.Ratio, Contact.Time..s., RSI..JH..Flight.Time..Contact.Time., RSI..Flight.Time.Contact.Time., Drop.Height..cm.)%>%
  mutate(RSI..m.per.s. = (Jump.Height..Flight.Time...cm./100)/Contact.Time..s.,
         Date = as.Date(Date, '%m/%d/%Y')) %>%
  group_by(Date, anon_id) %>%
  mutate(Weight..kg. = mean(na.omit(Weight..kg.)),
         Player.Average.RSI..m.per.s. = mean(RSI..m.per.s.)) %>%
  ungroup() %>%
  group_by(Date) %>%
  mutate(Team.Average.RSI..m.per.s. = mean(na.omit(RSI..m.per.s.)))

# Remove columns where all data points are NA
kinexon_session_clean <- kinexon_session |>
  select(where(~ !all(is.na(.))))


# Convert if Date column is NOT already in Date format
if (!inherits(kinexon_session_clean$Date, "Date")) {
  kinexon_session_clean <- kinexon_session_clean %>%
    mutate(Date = as.Date(Date, format = "%m/%d/%Y"))
}
# Confirm it's now in Date format
class(kinexon_session_clean$Date)

[1] "Date"

# Filter to only include data from this season
kinexon_session_clean <- kinexon_session_clean |>
  filter(Date >= as.Date("2024-6-01"))

##Analysis

Question 1

Are changes in RSI related to team game performance?

Team Wins and Losses

#filter ForceDecks data set into only the most recent season and only players that played in most of the games
current_season <- VALD_ForceDecks_clean %>%
  filter(Date >= as.Date("2024-09-20", '%Y-%m-%d')) %>%
  filter(!(anon_id %in% c("ID_10", "ID_14", "ID_15", "ID_18", "ID_54", "ID_64")))

#looking at the average RSI as well as the adjusted variance
#functions to calculate mean and variance of time series found in "Time-Series-Functions.Rmd"
team_mean <- mean(current_season$Team.Average.RSI..m.per.s.)
team_var <- as.numeric(sample_mean_and_variance(current_season$Team.Average.RSI..m.per.s.)[2])

#looking at the team average RSI over the course of the season, the middle horizontal line represents the mean and the gold lines represent one standard deviation away from the mean
ggplot(data = current_season ,aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line(linewidth = 0.75) +
  geom_hline(yintercept = team_mean) +
  geom_hline(yintercept = team_mean + sqrt(team_var), color="#CFB87C") +
  geom_hline(yintercept = team_mean - sqrt(team_var), color="#CFB87C")

#looking at the adjusted correlations of the team averages to see if there is any correlation with time
cor(x = current_season$X, y = current_season$Team.Average.RSI..m.per.s., method = "kendall")
cor(x = current_season$X, y = current_season$Team.Average.RSI..m.per.s., method = "spearman")

The plot above shows the team’s average RSI per week throughout the 2024-2025 season. This plot shows that there are large variations in team’s RSI values throughout the season but there does not seem to be any statistically significant underlying trend in the data. The lines in the plot represent the mean and one standard deviation from the mean throughout the season. As you can see, there is a lot of variability throughout the data. The patterns that appear though could potentially suggest an underlying ARIMA or other model for the data.

#calculating mean RSI values for each day and attaching them to correct dates
Team.means <- unique(current_season$Team.Average.RSI..m.per.s.)
Dates <- unique(current_season$Date)
IDs <- unique(current_season$anon_id)
Team.Change <- rep(0, 25)
current_season$Change.Team.Average.RSI..m.per.s. <- rep(NA, 580)

#calculating the difference in RSI from the previous measurement and attaching as a new column in the current season data set
for(i in 2:25){
  Team.Change[i] = Team.means[i] - Team.means[i+1]
  current_season[current_season$Date == Dates[i],19]<- Team.Change[i]
}

#making binary column of whether RSI increased or not
current_season <- current_season %>%
  mutate(Team.RSI.Increase = ifelse(Change.Team.Average.RSI..m.per.s. > 0, 1, 0))


#calculating the difference in RSI from the previous measurement for each player and attaching as a new column in the current season data set
current_season$Change.Player.Average.RSI..m.per.s. <- rep(0,580) #making new column

for(i in 1:11){ #cycle through players, calculate difference from last measurement
  Dates <- unique(current_season[current_season$anon_id == IDs[i],]$Date)
  Player.Means <- unique(current_season[current_season$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.) #taking measurements throughout season
  Player.Change <- rep(0,length(Player.Means))
  ifelse(i != 14, Player.Change <- rep(0,length(Player.Means)), 0) #don't need this line any more but I am not going to delete it

  for(j in 1:length(Player.Change)){ #calculates changes from last measurements
    Player.Change[j] <- Player.Means[j] - Player.Means[j+1]
    }
  #attaching RSI changes in new column  
  for(k in 1:length(Player.Change)){
    current_season[current_season$Date == Dates[k] & current_season$anon_id == IDs[i],]$Change.Player.Average.RSI..m.per.s. = Player.Change[k]
  }
}

#making indicator variable for whether RSI increased from last measurement or not
current_season <- current_season %>%
  mutate(Player.RSI.Increase = ifelse(Change.Player.Average.RSI..m.per.s. > 0, 1, 0))

Maybe Delete

#represents games throughout the season, will merge into one data set
Win <- c(1,1,1,1,1,0,1,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,1,1,1,0,0)
All_Dates <- as.Date(c("2024-10-19","2024-11-4", "2024-11-8", "2024-11-13", "2024-11-17", "2024-11-25","2024-11-26", "2024-11-27", "2024-12-2", "2024-12-7", "2024-12-13", "2024-12-21", "2024-12-30", "2025-1-4", "2025-1-8", "2025-1-12", "2025-1-15", "2025-1-18", "2025-1-21", "2025-1-25", "2025-1-28", "2025-2-2", "2025-2-5", "2025-2-8", "2025-2-11", "2025-2-15", "2025-2-18", "2025-2-22", "2025-2-24", "2025-3-2", "2025-3-5", "2025-3-8", "2025-3-11","2025-3-12", "2025-3-13", "2025-4-1"),"%Y-%m-%d")
Regular <- as.Date(c("2024-10-19","2024-11-4", "2024-11-8", "2024-11-13", "2024-11-17", "2024-12-2", "2024-12-7", "2024-12-13", "2024-12-21", "2024-12-30", "2025-1-4", "2025-1-8", "2025-1-12", "2025-1-15", "2025-1-18", "2025-1-21", "2025-1-25", "2025-1-28", "2025-2-2", "2025-2-5", "2025-2-8", "2025-2-11", "2025-2-15", "2025-2-18", "2025-2-22", "2025-2-24", "2025-3-2", "2025-3-5", "2025-3-8"),"%Y-%m-%d")
Tournament <- as.Date(c("2024-11-25","2024-11-26", "2024-11-27"),"%Y-%m-%d")
Post_Season <- as.Date(c("2025-3-11","2025-3-12", "2025-3-13", "2025-4-1"),"%Y-%m-%d")

#summarizes all games throughout the season
record <- data.frame(Date = All_Dates, Win = Win) %>%
  mutate(Type = case_when(
    Date %in% Regular ~ "Regular",
    Date %in% Tournament ~ "Tournament",
    Date %in% Post_Season ~ "Post-Season",
    TRUE ~ "NA"
  ))

#plotting regular game outcomes against RSI values
ggplot(current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line(linewidth = 0.75) +
  geom_vline(xintercept = record$Date[record$Type == "Regular"], linewidth = 0.2, color = ifelse(record$Win[record$Type == "Regular"] == 1, "palegreen3", "red3"))

#plotting tournaments or post season games with average team RSI values
ggplot(current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line(linewidth = 0.75) +
  geom_vline(xintercept = record$Date[record$Type %in% c("Tournament","Post-Season")], linewidth = 0.2, color = ifelse(record$Win[record$Type %in% c("Tournament", "Post-Season")] == 1, "palegreen3", "red3"))

These plots show that there may not be much predictive power for team game performance from team RSI values themselves. Instead, it may be useful to look at the changes in RSI before each game or look at a normalized plot of the team’s average RSI to get a look at what it would look like without as much noise.

The data was normalized with a moving average process which removed some of the noise present in the data. For this plot, the green lines represent a game in which the team’s stats were higher than the median for steals, assists, and rebounds. Red lines represent games where the total team sum of steals, assists, and rebounds are lower than the season median. Again, There does not seem to be a logical association between team average RSI and in game statistics. Instead, we find that games that had better in game performance are associated with lower RSI recordings before the game and games that had worse in game performance are associated with higher RSI recording just before the game. It’s important to note that there gaps between measurements that could have missed important spikes or dips in RSI for the team that would give us a more clear understanding of the relationship. Previous research that was able to identify significant relationships between in game performance and RSI measurements were able to collect it before and after each event.

#plotting the change in team average RSI throughout the season
ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0)


#plotting the changes in team average RSI throughout the season with green games being games where the team performed better than their median performance in explosiveness metrics and red being games where the team performed below median for explosiveness metrics
ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0) +
  geom_vline(xintercept = stats_clean$Date[stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1], color = ifelse(stats_clean[(stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1),]$Team.Good.Game == 1, "palegreen3","red3"))
#  geom_text(data = stats_clean, x = Date, y = 0.19, label = Date, angle = 90,size = 2.75)

This plot shows the change in team average RSI measurements throughout the season. Each data point represents how much the average RSI changed from the previous measurement. This plot helps to make sense of how the RSI of the team changed throughout the season. We can see that while increases were more frequent, they were smaller than the drops in RSI that occurred throughout the season. With this though, we can see that there is potentially a relationship between changes in RSI and in team game performance. This is because visually, it seems that games where the team had a “good” game tend to have an increase in RSI measured right before the game and games that are “bad” tend to have a decrease in RSI right before the game or a less dramatic increase in RSI than games that are “good”.

Maybe Delete

ggplot(data = current_season, aes(x=Date, y=Team.Average.RSI..m.per.s.)) +
  geom_line(aes(color = Team.RSI.Increase)) +
  geom_vline(xintercept = record$Date[record$Win == 1], linewidth = 0.2, color = "palegreen3")

ggplot(data = current_season, aes(x=Date, y=Team.Average.RSI..m.per.s.)) +
  geom_line(aes(color = Team.RSI.Increase), linewidth = 0.75) +
  geom_vline(xintercept = record$Date[record$Win == 0], linewidth = 0.2, color = "red3")

Team Explosiveness Metrics

ggplot(data = current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Points.Scored > median(stats_clean$Team.Points.Scored), "palegreen3", "red3"), linewidth = 0.2)


ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0) +
   geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Points.Scored > median(stats_clean$Team.Points.Scored), "palegreen3", "red3"), linewidth = 0.2)

This plot represents if the team scored higher or lower than the median team score in that game. Green represents higher and red represents lower than the median. This plot suggests that lower than median scores are associated with drops in RSI before the game (found in the middle of the season). As the season goes on, we can see that even though the team’s RSI goes up, the team has a consistently lower than median score. This could be indicative of trouble recovering from fatigue found in earlier research.

#Looking at games that had team total rebounds above or below median team season value
ggplot(data = current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Rebounds > median(stats_clean$Team.Rebounds), "palegreen3", "red3"), linewidth = 0.2)

ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0)  +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Rebounds > median(stats_clean$Team.Rebounds), "palegreen3", "red3"), linewidth = 0.2)

#Looking at games that had team total steals above or below median team season value
ggplot(data = current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Steals > median(stats_clean$Team.Steals), "palegreen3", "red3"), linewidth = 0.2)

  ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0)  +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Steals > median(stats_clean$Team.Steals), "palegreen3", "red3"), linewidth = 0.2)

#Looking at games that had team total assists above or below median team season value
ggplot(data = current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Assists > median(stats_clean$Team.Assists), "palegreen3", "red3"), linewidth = 0.2)

ggplot(data = current_season, aes(x=Date, y=Change.Team.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept = 0)  +
  geom_vline(xintercept = stats_clean$Date, color = ifelse(stats_clean$Team.Assists > median(stats_clean$Team.Assists), "palegreen3", "red3"), linewidth = 0.2)

ggplot(data = current_season, aes(Date, Team.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date[stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1], color = ifelse(stats_clean[(stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1),]$Team.Good.Game == 1, "palegreen3","red3"))

Applying regularization techniques to team data

Dates <- unique(current_season$Date)
#smooth <- exponential_smoothing(unique(current_season$Team.Average.RSI..m.per.s.),2,0)
ma <- moving_average(unique(current_season$Team.Average.RSI..m.per.s.), 1)
#x53 <- five_three_x(unique(current_season$Team.Average.RSI..m.per.s.))

#current_season$Team.Average.RSI.Smooth <- rep(NA,580)
current_season$Team.Average.RSI.Moving.Average <- rep(NA, 580)
#current_season$Team.Average.RSI.53X <- rep(NA, 580)

for(i in 1:25){
  #current_season[current_season$Date == Dates[i],21]<- smooth[i]
  current_season[current_season$Date == Dates[i],23]<- ma[i]
  #current_season[current_season$Date == Dates[i],23]<- x53[i]
}

ggplot(data=current_season, aes(Date, Team.Average.RSI.Moving.Average)) +
  geom_line() +
geom_vline(xintercept = stats_clean$Date[stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1], color = ifelse(stats_clean[(stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1),]$Team.Good.Game == 1, "palegreen3","red3"))

Assessing effectiveness of RSI as predictor of in game success

#getting dates that had games on them
Game.Days <- unique(stats_clean$Date)

#numbering which game in the season it was
Game.Number <- seq(30, 1, -1)

#making empty columns to be filled with game number
stats_clean$Game.Number <- rep(NA, 315)
current_season$Game.Number <- rep(NA, 580)

#going through games data set and labeling which game was which throughout the season
## first, second, third, ... 30th
for(i in 1:30){
  stats_clean[stats_clean$Date == Game.Days[i],]$Game.Number <- Game.Number[i]
}

#going through practices data set and labeling the closest game that the practice was before
## will be used for grouping games together later
for(i in 1:30){
  current_season[current_season$Date < Game.Days[i],]$Game.Number = 31-i
}

#making one data set with everything together to assess effectiveness of RSI as a predictor
## grouping by the game played or led up to and the player id
Q1_data <- left_join(x=current_season, y=stats_clean, by=c("anon_id", "Game.Number"))
Q1_data <- Q1_data %>%
  filter(Team.Good.Game %in% c(0,1))
Q1_data <- Q1_data[13:540,]
Q1_data <- Q1_data[-c(468,469,514,515,516),]

#set seed to make reproducible, splitting data into training and testing sets (75-25 split)
set.seed(123)
rows <- sample(1:nrow(Q1_data), nrow(Q1_data)*0.75, replace=FALSE)
Q1_train <- Q1_data[rows,]
Q1_test <- Q1_data[-rows,]

#model to predict whether game was a good game or not based on team average RSI increase leading up to that game
Q1_model_good_game <- glm(Team.Good.Game~Change.Team.Average.RSI..m.per.s., data=Q1_train, family="binomial")
summary(Q1_model_good_game)

#using the model above to make prediction on whether it was a good game or not and calculating CER (using 0.5 as threshold, may adjust later)
### Changed threshold to 0.43 based off of ROC curve done below
Q1_test$Team.Good.Pred <- ifelse(predict.glm(Q1_model_good_game, Q1_test, type="response") >0.46,1,0)

Q1_test %>%
  summarize(CER = mean(na.omit(Team.Good.Pred != Team.Good.Game)))

#building binomial model to predict whether it was a bad team game or not
Q1_model_bad_game <- glm(Team.Bad.Game~Change.Team.Average.RSI..m.per.s., data = Q1_train, family = "binomial")
summary(Q1_model_bad_game)

#making prediction based on the model and calculating the CER for the model on test data
Q1_test$Team.Bad.Pred <- ifelse(predict.glm(Q1_model_bad_game, Q1_test, type="response") >0.05,1,0)

Q1_test %>%
  summarize(CER = mean(na.omit(Team.Bad.Pred != Team.Bad.Game)))

Cross Validating the Team Models and Calculating Sensitivity and Specificity Rates w/ ROC

set.seed(123)
#model to be used for cross validation and ROC curve
model_good_team_cv <- glm(Team.Good.Game ~ Change.Team.Average.RSI..m.per.s., data=Q1_data, family="binomial")

#Making predictions for ROC curve
#$Team.Good.Game.Pred <- predict(model_good_team_cv, type="response")

#making ROC curve
#pred <- prediction(predictions = Q1_data$Team.Good.Game.Pred,
#                   labels = Q1_data$Team.Good.Game)
#perf <- performance(pred,measure='tnr',x.measure='tpr')
#plot(perf, colorize=TRUE)


#cost function for CER in cross validation
cost <- function(obs, pred){
  mean((pred <= 0.46) & obs==1 | (pred > 0.46) & obs==0)
}

#cross validating model
good_team_cv <- cv.glm(data=Q1_data,glmfit=model_good_team_cv,cost,K=5)

#extract average error
good_team_cv$delta[1]

Q1_data$Team.Good.Game.Pred <- ifelse(predict(model_good_team_cv, type="response")>0.46,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Team.Good.Game.Pred, y=Q1_data$Team.Good.Game)

set.seed(123)
#model to be used for cross validation and ROC curve
model_bad_team_cv <- glm(Team.Bad.Game ~ Change.Team.Average.RSI..m.per.s., data=Q1_data, family="binomial")

#Making predictions for ROC curve
#Q1_data$Team.Bad.Game.Pred <- predict(model_bad_team_cv, type="response")

#making ROC curve
#pred <- prediction(predictions = Q1_data$Team.Bad.Game.Pred,
#                   labels = Q1_data$Team.Bad.Game)
#perf <- performance(pred,measure='fnr',x.measure='fpr')
#plot(perf, colorize=TRUE)


#cost function for CER in cross validation
cost <- function(obs, pred){
  mean((pred <= 0.05) & obs==1 | (pred > 0.05) & obs==0)
}

#cross validating model
bad_team_cv <- cv.glm(data=Q1_data,glmfit=model_bad_team_cv,cost,K=5)

#extract average error
bad_team_cv$delta[1]

Q1_data$Team.Bad.Game.Pred <- ifelse(predict(model_bad_team_cv, type="response")>0.05,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Team.Bad.Game.Pred, y=Q1_data$Team.Bad.Game)

In order to limit variability and mimic past research, binary variables were created. Other research used top running speed in games to measure explosiveness of an athlete at a particular time and were compared to median top in game running speed throughout the season to define success. The data sets given did not have in game running speed so a combination of steals, assists, and rebounds were used as a proxy for a measure of explosiveness. A binary variable was created that measured if the team performed above or below median in all three metrics (1) or if one or more of the listed metrics were below the season team median (0). This is what will be predicted through changes in RSI of the team throughout the season. Shown above are the cross validated predicted error rates for the models that are used to predict whether the team will play above their median game play or below their median game play in a given game.

For the first model, it is trying to predict whether or not the team will have their steals, rebounds, and assists, all above their season median for a given game where the only predictor is the change in RSI measured the previous day. This model has a cross validated classification error rate of around 0.4 meaning that it will make a wrong prediction around 40% of the time. Out of the misclassifications, the model will generally tend to predict that a game will not be above median when it is. In fact, the model will predict that a game will not be above median when it actually is, about 3 times as often as it will correctly predict an above median game will be above median. Due to this major lack of sensitivity, it’s somewhat safe to say that changes in RSI may not be a good predictor for in game success based on the constraints made in this analysis.

The second model is trying to predict whether or not the team will have the total rebounds, assists, and steals for the game all below median or not. With this model, the cross validated classification error rate is around 0.2. This suggests that the model will make an incorrect classification around 20% of the time. Out of the incorrect classifications all of them were times in which in the model predicted for the team to perform below their median in explosiveness metrics in the game but the team performed above their median in at least one metric. This suggests that this model is making predictions that the team will perform “poorly” or below median more frequently than what is actually happening. In fact, when the model predicts that the team will perform below their season median game play, the predictions will be incorrect more often than correct. It’s important to note that this model correctly identified all of the games in which the team performed below median. This suggests that there is predictive power in changes in RSI for in game performance that is measured to be below median for the team.

For the data points that were misclassified by the model, it may be useful to look at the metrics for those games and identify how the team performed. Since this analysis defined a bad game as a game in which the team performed below median in steals, assists, and rebounds or not, it may be useful to see which metrics they performed below median in games that were predicted to be “bad games” but where not below median in all three metrics. The model could have potentially identified games that were poor but were not poor enough per say to be classified as a “bad game” in this analysis.

Question 2

Are changes in RSI related to individual statistical game performance?

#Looking at games that had above or below median rebounds throughout the season
ggplot(data = current_season[current_season$anon_id == "ID_42",], aes(Date, Player.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date[stats_clean$anon_id == "ID_42"], color = ifelse(stats_clean[stats_clean$anon_id == "ID_42",]$Rebounds >= median(stats_clean[stats_clean$anon_id == "ID_42",]$Rebounds), "palegreen3", "red3"))

#Looking at games that had above or below median steals throughout the season
ggplot(data = current_season[current_season$anon_id == "ID_42",], aes(Date, Player.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date[stats_clean$anon_id == "ID_42"], color = ifelse(stats_clean[stats_clean$anon_id == "ID_42",]$Steals >= median(stats_clean[stats_clean$anon_id == "ID_42",]$Steals), "palegreen3", "red3"))

#Looking at games that had above or below median assists throughout the season
ggplot(data = current_season[current_season$anon_id == "ID_42",], aes(Date, Player.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date[stats_clean$anon_id == "ID_42"], color = ifelse(stats_clean[stats_clean$anon_id == "ID_42",]$Assists >= median(stats_clean[stats_clean$anon_id == "ID_42",]$Assists), "palegreen3", "red3"))

ggplot(data = current_season[current_season$anon_id == "ID_42",], aes(Date, Player.Average.RSI..m.per.s.)) +
  geom_line() +
  geom_vline(xintercept = stats_clean$Date[stats_clean$anon_id == "ID_42" & (stats_clean$Player.Good.Game == 1 | stats_clean$Player.Bad.Game ==1)], color = ifelse(stats_clean[stats_clean$anon_id == "ID_42" & (stats_clean$Player.Good.Game == 1 | stats_clean$Player.Bad.Game ==1),]$Player.Good.Game == 1, "palegreen3","red3"))

ggplot(data = current_season[current_season$anon_id == "ID_42",], aes(Date, Change.Player.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_hline(yintercept=0) +
  geom_vline(xintercept = stats_clean$Date[stats_clean$anon_id == "ID_42" & (stats_clean$Player.Good.Game == 1 | stats_clean$Player.Bad.Game ==1)], color = ifelse(stats_clean[stats_clean$anon_id == "ID_42" & (stats_clean$Player.Good.Game == 1 | stats_clean$Player.Bad.Game ==1),]$Player.Good.Game == 1, "palegreen3","red3"))

Applying regularization techniques to player data

smooth <- exponential_smoothing(unique(current_season$Team.Average.RSI..m.per.s.),2,0)
ma <- moving_average(unique(current_season$Team.Average.RSI..m.per.s.), 1)
x53 <- five_three_x(unique(current_season$Team.Average.RSI..m.per.s.))

current_season$Team.Average.RSI.Smooth <- rep(NA,580)
current_season$Team.Average.RSI.Moving.Average <- rep(NA, 580)
current_season$Team.Average.RSI.53X <- rep(NA, 580)

for(i in 1:25){
  current_season[current_season$Date == Dates[i],21]<- smooth[i]
  current_season[current_season$Date == Dates[i],22]<- ma[i]
  current_season[current_season$Date == Dates[i],23]<- x53[i]
}

ggplot(data=current_season, aes(Date, Team.Average.RSI.Moving.Average)) +
  geom_line() +
geom_vline(xintercept = stats_clean$Date[stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1], color = ifelse(stats_clean[(stats_clean$Team.Good.Game ==1 | stats_clean$Team.Bad.Game ==1),]$Team.Good.Game == 1, "palegreen3","red3"))

Assessing effectiveness of RSI as predictor of in game success

Good Games Models

#building completely pooled model for player RSI changes and good games
Q2_model_good_game_pooled <- glm(Player.Good.Game~Change.Player.Average.RSI..m.per.s., data=Q1_train, family = "binomial")
summary(Q2_model_good_game_pooled)

#making predictions with the pooled model for good games
Q1_test$Player.Good.Pooled.Pred <- ifelse(predict(Q2_model_good_game_pooled, Q1_test)>0.5,1,0)

#calculating the CER of the model on the test data
#Q1_test %>%
# summarize(CER = mean(Player.Good.Game != Player.Good.Pooled.Pred))

#making unpooled model for players
Q2_model_good_game_unpooled <-glmer(Player.Good.Game~Change.Player.Average.RSI..m.per.s. + (Change.Player.Average.RSI..m.per.s.|anon_id), data=Q1_train, family = "binomial")
summary(Q2_model_good_game_unpooled)

#making predictions for unpooled model on the test data (calculated with probability)
Q1_test$Player.Good.Unpooled.Pred <- ifelse(predict(Q2_model_good_game_unpooled, Q1_test, type="response")>0.46,1,0)

#calculating CER of unpooled model
Q1_test %>%
  summarize(CER = mean(Player.Good.Unpooled.Pred != Player.Good.Game))

#building completely pooled model for player RSI changes and bad games
Q2_model_bad_game_pooled <- glm(Player.Bad.Game~Change.Player.Average.RSI..m.per.s., data=Q1_train, family = "binomial")
summary(Q2_model_bad_game_pooled)

#making predictions for bad games with pooled model on test data
Q1_test$Player.Bad.Pooled.Pred <- ifelse(predict(Q2_model_bad_game_pooled, Q1_test)>0.5,1,0)

#calculating CER for bad games with pooled model
#Q1_test %>%
#  summarize(CER = mean(Player.Bad.Game != Player.Bad.Pooled.Pred))

#making unpooled model for bad games
Q2_model_bad_game_unpooled <-glmer(Player.Bad.Game~Change.Player.Average.RSI..m.per.s. + (Change.Player.Average.RSI..m.per.s.|anon_id), data=Q1_train, family = "binomial")
summary(Q2_model_bad_game_unpooled)

Q1_test$Player.Bad.Unpooled.Pred <- ifelse(predict(Q2_model_bad_game_unpooled, Q1_test, type="response")>0.03,1,0)

Q1_test %>%
  summarize(CER = mean(Player.Bad.Unpooled.Pred != Player.Bad.Game))

Cross Validating Errors for Good and Bad Games with ROC Curves

set.seed(123)

#model to be used for cross validation
model_good_player_pooled_cv <- glm(Player.Good.Game ~ Change.Player.Average.RSI..m.per.s., data=Q1_data, family="binomial")


#Q1_data$Player.Good.Game.Pooled.Pred <- predict(model_good_player_pooled_cv, type="response")

#pred <- prediction(predictions = Q1_data$Player.Good.Game.Pooled.Pred,
#                   labels = Q1_data$Player.Good.Game)
#perf <- performance(pred,measure='tnr',x.measure='tpr')
#plot(perf, colorize=TRUE)
#abline(a=0, b=1)

cost <- function(obs, pred){
  mean((pred <=.49) &obs==1 | (pred > 0.49) & obs==0)
}

#cross validating model
good_player_cv <- cv.glm(data=Q1_data,glmfit=model_good_player_pooled_cv,cost,K=5)

#extract average error
good_player_cv$delta[1]

Q1_data$Player.Good.Game.Pooled.Pred <- ifelse(predict(model_good_player_pooled_cv, type="response")>0.49,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Player.Good.Game.Pooled.Pred, y=Q1_data$Player.Good.Game)

set.seed(123)

#model to be used for cross validation
model_good_player_unpooled_cv <- glmer(Player.Good.Game~Change.Player.Average.RSI..m.per.s. + (Change.Player.Average.RSI..m.per.s.|anon_id), data=Q1_data, family = "binomial")


#Q1_data$Player.Good.Game.Unpooled.Pred <- predict(model_good_player_unpooled_cv, type="response")

#pred <- prediction(predictions = Q1_data$Player.Good.Game.Unpooled.Pred,
#                   labels = Q1_data$Player.Good.Game)
#perf <- performance(pred,measure='tnr',x.measure='tpr')
#plot(perf, colorize=TRUE)
#abline(a=0, b=1)


#making cost function for cross validation for unpooled model
cost <- function(y, yhat){
  mean(((yhat <= 0.46) & y==1) | ((yhat > 0.46) & y==0))
}


#cross validating the unpooled model
summary(cv(model_good_player_unpooled_cv, k=5, clusterVariables = "anon_id", seed=123, criterion=cost))


Q1_data$Player.Good.Game.Unpooled.Pred <- ifelse(predict(model_good_player_unpooled_cv, type="response")>0.46,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Player.Good.Game.Unpooled.Pred, y=Q1_data$Player.Good.Game)

set.seed(123)

#model to be used for cross validation
model_bad_player_unpooled_cv <- glmer(Player.Bad.Game~Change.Player.Average.RSI..m.per.s. + (1|anon_id) + (0+Change.Player.Average.RSI..m.per.s.|anon_id), data=Q1_data, family = "binomial")

#Q1_data$Player.Bad.Game.Unpooled.Pred <- predict(model_bad_player_unpooled_cv, type="response")

#pred <- prediction(predictions = Q1_data$Player.Bad.Game.Unpooled.Pred,
#                   labels = Q1_data$Player.Bad.Game)
#perf <- performance(pred,measure='tnr',x.measure='tpr')
#plot(perf, colorize=TRUE)
#abline(a=0, b=1)


#making cost function for cross validation for unpooled model
cost <- function(y, yhat){
  mean((yhat <=.05) & y==1 | (yhat > 0.05) & y==0)
}

summary(cv(model_bad_player_unpooled_cv, k=5, clusterVariables = "anon_id", seed=123, criterion=cost))

Q1_data$Player.Bad.Game.Unpooled.Pred <- ifelse(predict(model_bad_player_unpooled_cv, type="response")>0.05,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Player.Bad.Game.Unpooled.Pred, y=Q1_data$Player.Bad.Game)

set.seed(123)

#model to be used for cross validation
model_bad_player_cv <- glm(Player.Bad.Game ~ Change.Player.Average.RSI..m.per.s., data=Q1_data, family="binomial")

#Q1_data$Player.Bad.Game.Pooled.Pred <- predict(model_bad_player_cv, type="response")

#pred <- prediction(predictions = Q1_data$Player.Bad.Game.Pooled.Pred,
#                   labels = Q1_data$Player.Bad.Game)
#perf <- performance(pred,measure='tnr',x.measure='tpr')
#plot(perf, colorize=TRUE)
#abline(a=0, b=1)

cost <- function(obs, pred){
  mean((pred <=.03) & obs==1 | (pred > 0.03) & obs==0)
}

#cross validating model
bad_player_cv <- cv.glm(data=Q1_data,glmfit=model_bad_player_cv,cost,K=5)

#extract average error
bad_team_cv$delta[1]

Q1_data$Player.Bad.Game.Pooled.Pred <- if_else(predict(model_bad_player_cv, type="response")>0.03,1,0)

#looking at sensitivity and specificity
table(x=Q1_data$Player.Bad.Game.Pooled.Pred, y=Q1_data$Player.Bad.Game)

errors <- data.frame(Model = c("Good-Games-Pooled", "Good-Games-Unpooled", 
                                    "Bad-Games-Pooled", "Bad-Games-Unpooled"),
                     CER = c(0.4503817, 0.3740458   , 0.03053435, 0.2824427),
                     CER.CV = c(0.4837476, 0.456979, 0.2045889, 0.03250478),
                     Classification.Threshold = c(0.49, 0.46, 0.03, 0.05))

kable(errors)

Based on the errors that come from the cross validation of all of the models, we can see that the unpooled models outperform the pooled models. By fitting a unique model to each player we were able to get better predictions for in game success than fitting a common model to all of the players. This suggests that RSI may have differing predictive power for in game performance for each player. This means that changes in RSI may be a more powerful predictor for in game success than for other athletes on the team. This also suggests that some information may be lost when averaging changes in RSI for the team. One thing to note about these models is that the metrics shown above are predicted rates in which the model will make a wrong prediction. It’s important to know that each model gives the predicted probability that the given observation is a “good game” for the first two models or the probability that a given observation is a “bad game” for the last two models. Normally the threshold to consider a prediction to be a 1 or a 0 is if the probability is above or below 0.5 (rounded to the nearest whole number). These models were given different classification thresholds. This was for computational reasons since “bad games” for players are very rare and will therefore be harder to identify. For the models that are predicting whether a game will be a “good game” or not, they don’t have a lot of predictive power.

Question 3:

Is the previous week’s load related to RSI?

Recalculating Gameday Minus Columns

# Recalculating the Gameday.Minus.N.Acceleration.Load Columns

# Create a df with relevant columns and add column specifying if that day was a game or not
daily_load <- kinexon_session_clean |>
  mutate(Game.Day = if_else(Type == "Game", 1, 0)) |>
  filter(!is.na(Daily.Total.Accel.Load.Accum)) |>
  select(anon_id, Date, Daily.Total.Accel.Load.Accum, Game.Day) |>
  arrange(anon_id, Date)

# Get all game day dates per player as a lookup table
game_days <- daily_load |>
  filter(Game.Day == 1) |>
  select(anon_id, Game.Date = Date)
    

# Add column Day.Type for GD-1 - GD-4 and gamedays
daily_load <- daily_load |>
  left_join(game_days, by = "anon_id", relationship = "many-to-many") |>
  mutate(
    Days.Until.Game = as.numeric(difftime(Game.Date, Date, units = "days")),
    Day.Type = case_when(
      Days.Until.Game == 1 ~ "GD-1",
      Days.Until.Game == 2 ~ "GD-2",
      Days.Until.Game == 3 ~ "GD-3",
      Days.Until.Game == 4 ~ "GD-4",
      Game.Day == 1 ~ "GameDay",
      TRUE ~ NA_character_
    )
  ) |>
  filter(!is.na(Day.Type)) |>    # keep only rows within 4 days before game
  select(-Game.Date, -Days.Until.Game)

# Add columns for Gameday.Minus.N.Acceleration.Load
Gameday_Minus <- daily_load |>
  mutate(
    Gameday.Minus.1.Acceleration.Load = if_else(Day.Type == "GD-1", Daily.Total.Accel.Load.Accum, NA_real_),
    Gameday.Minus.2.Acceleration.Load = if_else(Day.Type == "GD-2", Daily.Total.Accel.Load.Accum, NA_real_),
    Gameday.Minus.3.Acceleration.Load = if_else(Day.Type == "GD-3", Daily.Total.Accel.Load.Accum, NA_real_),
    Gameday.Minus.4.Acceleration.Load = if_else(Day.Type == "GD-4", Daily.Total.Accel.Load.Accum, NA_real_),
    GameDay.Acceleration.Load = if_else(Day.Type == "GameDay", Daily.Total.Accel.Load.Accum, NA_real_)
  )

# Delete duplicate rows
Gameday_Minus <- Gameday_Minus |>
  distinct()

# Create columns for Average.Gameday.Minus.N.Acceleration.Load

# Calculate and add columns for running averages
Gameday_Minus <- Gameday_Minus |>
  arrange(anon_id, Date) |>
  group_by(anon_id) |>
  mutate(
    # Cumulative count and sum for GD-1
    cum_count_gd1 = cumsum(!is.na(Gameday.Minus.1.Acceleration.Load)),
    cum_sum_gd1 = cumsum(replace_na(Gameday.Minus.1.Acceleration.Load, 0)),
    Average.Gameday.Minus.1.Acceleration.Load = if_else(
      cum_count_gd1 > 0,
      cum_sum_gd1 / cum_count_gd1,
      NA_real_
    ),
    
    # Cumulative count and sum for GD-2
    cum_count_gd2 = cumsum(!is.na(Gameday.Minus.2.Acceleration.Load)),
    cum_sum_gd2 = cumsum(replace_na(Gameday.Minus.2.Acceleration.Load, 0)),
    Average.Gameday.Minus.2.Acceleration.Load = if_else(
      cum_count_gd2 > 0,
      cum_sum_gd2 / cum_count_gd2,
      NA_real_
    ),
    
    # Cumulative count and sum for GD-3
    cum_count_gd3 = cumsum(!is.na(Gameday.Minus.3.Acceleration.Load)),
    cum_sum_gd3 = cumsum(replace_na(Gameday.Minus.3.Acceleration.Load, 0)),
    Average.Gameday.Minus.3.Acceleration.Load = if_else(
      cum_count_gd3 > 0,
      cum_sum_gd3 / cum_count_gd3,
      NA_real_
    ),
    
    # Cumulative count and sum for GD-4
    cum_count_gd4 = cumsum(!is.na(Gameday.Minus.4.Acceleration.Load)),
    cum_sum_gd4 = cumsum(replace_na(Gameday.Minus.4.Acceleration.Load, 0)),
    Average.Gameday.Minus.4.Acceleration.Load = if_else(
      cum_count_gd4 > 0,
      cum_sum_gd4 / cum_count_gd4,
      NA_real_
    )
  ) |>
  ungroup() |>
  select(-starts_with("cum_count_"), -starts_with("cum_sum_")) # clean up columns

# Calculate running averages of Gameday.Acceleration.Load and add column Average.Gameday.Acceleration.Load

Gameday_Minus <- Gameday_Minus |>
  arrange(anon_id, Date) |>
  group_by(anon_id) |>
  mutate(
    cum_count_gd = cumsum(!is.na(GameDay.Acceleration.Load)),
    cum_sum_gd = cumsum(replace_na(GameDay.Acceleration.Load, 0)),
    Average.GameDay.Acceleration.Load = if_else(
      cum_count_gd > 0,
      cum_sum_gd / cum_count_gd,
      NA_real_
    )
  ) |>
  ungroup() |>
  select(-cum_count_gd, -cum_sum_gd) # clean up columns

# Filter to only have data from this season
RSI <- VALD_ForceDecks_clean |>
  filter(Date >= as.Date("2024-10-10")) |>
  select(anon_id, Date, RSI..m.per.s., Trial) |>
  arrange(anon_id, Date) |>
  group_by(anon_id, Date) |>
  mutate(Daily.Avg.RSI..m.per.s = mean(RSI..m.per.s.)) |>
  ungroup()


# Plot RSI over season
ggplot(RSI, aes(x=Date, y=Daily.Avg.RSI..m.per.s, color=anon_id)) +
  geom_point() +
  geom_smooth(method="lm", se=FALSE) +
  labs(title = "Athlete RSI Scores Throughout Season")

`geom_smooth()` using formula = 'y ~ x'

# Create load data frame containing columns relevant to load
load <- kinexon_session_clean |>
  filter(Date >= as.Date("2024-10-10")) |>
  select(anon_id, Date, Daily.Total.Accel.Load.Accum, Jump.Load) |>
  filter(!is.na(Daily.Total.Accel.Load.Accum))


# Look into Daily Acceleration Load for one athlete
id_40_load <- load |>
  filter(anon_id == "ID_40")
# Plot of ID_40's Daily Accel Load
ggplot(data = id_40_load, aes(x = Date, y = Daily.Total.Accel.Load.Accum)) +
  geom_point() +
  geom_line() +
  labs(title = "ID_40 Daily Acceleration Load")

# Plot of all athlete's daily accel load throughout season
ggplot(data = load, aes(x=Date, y=Daily.Total.Accel.Load.Accum, color=anon_id)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "2024-25 Baskteball Athlete's Daily Acceleration Load")

`geom_smooth()` using formula = 'y ~ x'

# Plot of all athlete's daily jump load throughout season
ggplot(data = load, aes(x=Date, y=Jump.Load, color=anon_id)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  labs(title = "2024-25 Baskteball Athlete's Daily Jump Load")

`geom_smooth()` using formula = 'y ~ x'

The graphs seem to show that for most athletes, RSI scores increased throughout the season, while daily acceleration load decreased and jump load about stayed the same.

# Calculate the average and total acceleration and jump load in the last seven days.
# Last 7 days not last 7 sessions
# Include that day in the last 7.
rolling_7day <- load |>
  filter(!is.na(Daily.Total.Accel.Load.Accum)) |>
  arrange(anon_id, Date) |>
  group_by(anon_id) |>
  group_split() |>
  map_dfr(~{
    athlete_data <- .
    map_dfr(seq_len(nrow(athlete_data)), function(i) {
      current_row <- athlete_data[i, ]
      past_week <- athlete_data |>
        filter(Date <= current_row$Date & Date > current_row$Date - days(7)) # (<=) includes that day in the last 7.
      
      current_row |>
        mutate(
          Last.Week..Avg..Accel.Load = mean(past_week$Daily.Total.Accel.Load.Accum, na.rm = TRUE),
          Last.Week..Total..Accel.Load = sum(past_week$Daily.Total.Accel.Load.Accum, na.rm = TRUE),
          Last.Week..Avg..Jump.Load = mean(past_week$Jump.Load, na.rm = TRUE),
          Last.Week..Total..Jump.Load = sum(past_week$Jump.Load, na.rm = TRUE),
          Days.Used = nrow(past_week)
        )
    })
  })

RSI_daily_avg <- RSI |>
  distinct(anon_id, Date, Daily.Avg.RSI..m.per.s)

# Merge the RSI and load data frames
load_and_RSI <- RSI_daily_avg |>
  left_join(
    rolling_7day,
    by = "anon_id",
    relationship = "many-to-many"
  ) |>
  filter(Date.y <= Date.x & Date.y >= Date.x - days(7)) |> # <= includes same day. On days when RSI and accel load are measured, accel load is measured first, meaning we should include it.
  group_by(anon_id, Date.x) |>
  slice_max(Date.y, n = 1) |>  # Keep the most recent load data before RSI
  ungroup() |>
  rename(
    RSI_Date = Date.x,
    Load_Date = Date.y,
    RSI = Daily.Avg.RSI..m.per.s
  )

# Average Accel Load:
# Pearson:
cor(load_and_RSI$RSI, load_and_RSI$Last.Week..Avg..Accel.Load, use = "complete.obs")

[1] 0.007807349

# Spearman
cor(load_and_RSI$RSI, load_and_RSI$Last.Week..Avg..Accel.Load, method = "spearman", use = "complete.obs")

[1] 0.02250988

# Total Accel Load:
# Pearson:
cor(load_and_RSI$RSI, load_and_RSI$Last.Week..Total..Accel.Load, use = "complete.obs")

[1] 0.0424607

# Spearman
cor(load_and_RSI$RSI, load_and_RSI$Last.Week..Total..Accel.Load, method = "spearman", use = "complete.obs")

[1] 0.05468703

# Average Jump Load:
# Pearson:
cor(load_and_RSI$RSI, load_and_RSI$Last.Week..Avg..Jump.Load, use = "complete.obs")

[1] 0.2584746

# Spearman
cor(load_and_RSI$RSI, load_and_RSI$Last.Week..Avg..Jump.Load, method = "spearman", use = "complete.obs")

[1] 0.2313584

# Total Jump Load:
# Pearson:
cor(load_and_RSI$RSI, load_and_RSI$Last.Week..Total..Jump.Load, use = "complete.obs")

[1] 0.2377019

# Spearman
cor(load_and_RSI$RSI, load_and_RSI$Last.Week..Total..Jump.Load, method = "spearman", use = "complete.obs")

[1] 0.2298689

# At first glance, it seems like jump load correlates more with RSI than acceleration load

Jump Load seems to be have higher correlation to RSI. Lowest correlation are the Last.Week..Avg..Accel.Load and RSI.

# Last week's load vs. RSI
ggplot(load_and_RSI, aes(x = Last.Week..Avg..Accel.Load, y = RSI, color=anon_id)) +
  geom_point() +
  geom_smooth(method = "lm", se=FALSE) +
  labs(title = "Relationship Between Previous Week's Accel Load and RSI")

`geom_smooth()` using formula = 'y ~ x'

# Plot for each athlete
unique_ids <- unique(load_and_RSI$anon_id)

for (id in unique_ids) {
  p <- ggplot(filter(load_and_RSI, anon_id == id),
              aes(x = Last.Week..Total..Accel.Load, y = RSI)) +
    geom_point() +
    geom_smooth(method = "lm") +
    labs(title = paste("RSI vs. Prior Week's Accel Load:", id),
         x = "Total Load Last 7 Days",
         y = "RSI")
   print(p)
   
}

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

# Per athlete correlations

# ACCEL LOAD:
# Average 7 day accel load
per_athlete_corrs_avg <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, Last.Week..Avg..Accel.Load, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_avg, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Correlation Between Avg 7-Day Accel Load and RSI",
       x = "Athlete",
       y = "Correlation (r)")

# Total 7 day accel load
per_athlete_corrs_total <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, Last.Week..Total..Accel.Load, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_total, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Correlation Between Total 7-Day Accel Load and RSI",
       x = "Athlete",
       y = "Correlation (r)")

# JUMP LOAD
# Average 7 day jump load
per_athlete_corrs_avg <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, Last.Week..Avg..Jump.Load, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_avg, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Correlation Between Avg 7-Day Jump Load and RSI",
       x = "Athlete",
       y = "Correlation (r)")

# Total 7 day jump load
per_athlete_corrs_total <- load_and_RSI |>
  group_by(anon_id) |>
  summarize(
    r = cor(RSI, Last.Week..Total..Jump.Load, use = "complete.obs"),
    n = n()
  ) |>
  arrange(r)

ggplot(per_athlete_corrs_total, aes(x = reorder(anon_id, r), y = r)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Correlation Between Total 7-Day Jump Load and RSI",
       x = "Athlete",
       y = "Correlation (r)")

# This allows us to see which athlete's RSI scores are most effected by the prior week's load

# Scatter plot of Total Jump Load and RSI
ggplot(load_and_RSI, aes(Last.Week..Total..Jump.Load, RSI)) +
  geom_point() +
  geom_smooth() +
  labs(title = "RSI vs. Total Weekly Jump Load")

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# Scatter plot of Avg Jump Load and RSI
ggplot(load_and_RSI, aes(Last.Week..Avg..Jump.Load, RSI)) +
  geom_point() +
  geom_smooth() +
  labs(title = "RSI vs. Avg Weekly Jump Load")

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

# Scatter plot of Total Jump Load and RSI with athletes
ggplot(load_and_RSI, aes(Last.Week..Total..Jump.Load, RSI, color=anon_id)) +
  geom_point() +
  geom_smooth(method="lm", se=FALSE) +
  labs(title = "RSI vs. Total Weekly Jump Load Per Athlete")

`geom_smooth()` using formula = 'y ~ x'

# Scatter plot of Avg Jump Load and RSI with athletes
ggplot(load_and_RSI, aes(Last.Week..Avg..Jump.Load, RSI, color=anon_id)) +
  geom_point() +
  geom_smooth(method="lm", se=FALSE) +
  labs(title = "RSI vs. Avg Weekly Jump Load Per Athlete")

`geom_smooth()` using formula = 'y ~ x'

set.seed(42)

# Create grouped 5-fold CV splits by anon_id
cv_splits <- group_vfold_cv(load_and_RSI, group = anon_id, v = 5)

# Mixed Effects Model with standardized predictors
cv_mixed <- function(split) {
  train <- training(split)
  test <- testing(split)

  # Manually standardize using training stats
  mean_accel <- mean(train$Last.Week..Avg..Accel.Load, na.rm = TRUE)
  sd_accel   <- sd(train$Last.Week..Avg..Accel.Load, na.rm = TRUE)

  mean_jump  <- mean(train$Last.Week..Avg..Jump.Load, na.rm = TRUE)
  sd_jump    <- sd(train$Last.Week..Avg..Jump.Load, na.rm = TRUE)

  train <- train %>%
    mutate(
      z_avg_accel = (Last.Week..Avg..Accel.Load - mean_accel) / sd_accel,
      z_avg_jump  = (Last.Week..Avg..Jump.Load - mean_jump) / sd_jump
    )

  test <- test %>%
    mutate(
      z_avg_accel = (Last.Week..Avg..Accel.Load - mean_accel) / sd_accel,
      z_avg_jump  = (Last.Week..Avg..Jump.Load - mean_jump) / sd_jump
    )

  model <- lmer(RSI ~ z_avg_accel + z_avg_jump + (1 | anon_id), data = train)
  preds <- predict(model, newdata = test, allow.new.levels = TRUE)

  rmse_vec(truth = test$RSI, estimate = preds)
}

# GAM Model with standardized predictors
cv_gam <- function(split) {
  train <- training(split)
  test <- testing(split)

  # Manually standardize using training stats
  mean_accel <- mean(train$Last.Week..Avg..Accel.Load, na.rm = TRUE)
  sd_accel   <- sd(train$Last.Week..Avg..Accel.Load, na.rm = TRUE)

  mean_jump  <- mean(train$Last.Week..Avg..Jump.Load, na.rm = TRUE)
  sd_jump    <- sd(train$Last.Week..Avg..Jump.Load, na.rm = TRUE)

  train <- train %>%
    mutate(
      z_avg_accel = (Last.Week..Avg..Accel.Load - mean_accel) / sd_accel,
      z_avg_jump  = (Last.Week..Avg..Jump.Load - mean_jump) / sd_jump
    )

  test <- test %>%
    mutate(
      z_avg_accel = (Last.Week..Avg..Accel.Load - mean_accel) / sd_accel,
      z_avg_jump  = (Last.Week..Avg..Jump.Load - mean_jump) / sd_jump
    )

  model <- gam(RSI ~ s(z_avg_accel) + s(z_avg_jump), data = train)
  preds <- predict(model, newdata = test)

  rmse_vec(truth = test$RSI, estimate = preds)
}

# Run CV for both models
mixed_rmses <- map_dbl(cv_splits$splits, cv_mixed)
gam_rmses   <- map_dbl(cv_splits$splits, cv_gam)

Error in `map_dbl()`:
ℹ In index: 1.
Caused by error in `gam()`:
! could not find function "gam"
Run `]8;;x-r-run:rlang::last_trace()rlang::last_trace()]8;;` to see where the error occurred.

# Fit fixed effects model on full data
avg_fixed_model <- lm(RSI ~ Last.Week..Avg..Accel.Load + Last.Week..Avg..Jump.Load + factor(anon_id),
                  data = load_and_RSI)

# View summary
summary(avg_fixed_model)

Call:
lm(formula = RSI ~ Last.Week..Avg..Accel.Load + Last.Week..Avg..Jump.Load + 
    factor(anon_id), data = load_and_RSI)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.64289 -0.11695  0.00447  0.13812  0.95503 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 2.374e+00  2.028e-01  11.702  < 2e-16 ***
Last.Week..Avg..Accel.Load -9.829e-04  2.655e-04  -3.702 0.000284 ***
Last.Week..Avg..Jump.Load   8.335e-06  4.074e-06   2.046 0.042179 *  
factor(anon_id)ID_14       -6.426e-01  1.744e-01  -3.684 0.000303 ***
factor(anon_id)ID_24        1.344e-02  1.554e-01   0.087 0.931143    
factor(anon_id)ID_40       -5.309e-01  1.531e-01  -3.467 0.000658 ***
factor(anon_id)ID_42        9.733e-02  1.558e-01   0.625 0.532970    
factor(anon_id)ID_5        -9.307e-01  1.528e-01  -6.090  6.6e-09 ***
factor(anon_id)ID_52       -5.781e-01  1.576e-01  -3.668 0.000322 ***
factor(anon_id)ID_6        -5.617e-01  1.554e-01  -3.614 0.000391 ***
factor(anon_id)ID_61       -1.253e-01  1.548e-01  -0.809 0.419570    
factor(anon_id)ID_62       -7.016e-02  1.614e-01  -0.435 0.664247    
factor(anon_id)ID_63       -5.128e-01  1.542e-01  -3.326 0.001066 ** 
factor(anon_id)ID_64       -1.338e-01  1.766e-01  -0.758 0.449681    
factor(anon_id)ID_66       -4.717e-01  1.494e-01  -3.157 0.001870 ** 
factor(anon_id)ID_74        2.652e-01  1.609e-01   1.648 0.101123    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2336 on 181 degrees of freedom
Multiple R-squared:  0.7177,    Adjusted R-squared:  0.6943 
F-statistic: 30.68 on 15 and 181 DF,  p-value: < 2.2e-16

set.seed(42)

# Create grouped 5-fold CV splits by anon_id
cv_splits <- group_vfold_cv(load_and_RSI, group = anon_id, v = 5)

# Mixed Effects Model with standardized total load predictors
cv_mixed_total <- function(split) {
  train <- training(split)
  test <- testing(split)

  # Manually standardize using training stats
  mean_accel <- mean(train$Last.Week..Total..Accel.Load, na.rm = TRUE)
  sd_accel   <- sd(train$Last.Week..Total..Accel.Load, na.rm = TRUE)

  mean_jump  <- mean(train$Last.Week..Total..Jump.Load, na.rm = TRUE)
  sd_jump    <- sd(train$Last.Week..Total..Jump.Load, na.rm = TRUE)

  train <- train %>%
    mutate(
      z_total_accel = (Last.Week..Total..Accel.Load - mean_accel) / sd_accel,
      z_total_jump  = (Last.Week..Total..Jump.Load - mean_jump) / sd_jump
    )

  test <- test %>%
    mutate(
      z_total_accel = (Last.Week..Total..Accel.Load - mean_accel) / sd_accel,
      z_total_jump  = (Last.Week..Total..Jump.Load - mean_jump) / sd_jump
    )

  model <- lmer(RSI ~ z_total_accel + z_total_jump + (1 | anon_id), data = train)
  preds <- predict(model, newdata = test, allow.new.levels = TRUE)

  rmse_vec(truth = test$RSI, estimate = preds)
}

# GAM Model with standardized total load predictors
cv_gam_total <- function(split) {
  train <- training(split)
  test <- testing(split)

  # Manually standardize using training stats
  mean_accel <- mean(train$Last.Week..Total..Accel.Load, na.rm = TRUE)
  sd_accel   <- sd(train$Last.Week..Total..Accel.Load, na.rm = TRUE)

  mean_jump  <- mean(train$Last.Week..Total..Jump.Load, na.rm = TRUE)
  sd_jump    <- sd(train$Last.Week..Total..Jump.Load, na.rm = TRUE)

  train <- train %>%
    mutate(
      z_total_accel = (Last.Week..Total..Accel.Load - mean_accel) / sd_accel,
      z_total_jump  = (Last.Week..Total..Jump.Load - mean_jump) / sd_jump
    )

  test <- test %>%
    mutate(
      z_total_accel = (Last.Week..Total..Accel.Load - mean_accel) / sd_accel,
      z_total_jump  = (Last.Week..Total..Jump.Load - mean_jump) / sd_jump
    )

  model <- gam(RSI ~ s(z_total_accel) + s(z_total_jump), data = train)
  preds <- predict(model, newdata = test)

  rmse_vec(truth = test$RSI, estimate = preds)
}

# Run CV for both models with total loads
mixed_rmses_total <- map_dbl(cv_splits$splits, cv_mixed_total)
gam_rmses_total   <- map_dbl(cv_splits$splits, cv_gam_total)

Error in `map_dbl()`:
ℹ In index: 1.
Caused by error in `gam()`:
! could not find function "gam"
Run `]8;;x-r-run:rlang::last_trace()rlang::last_trace()]8;;` to see where the error occurred.

Question 4

What is each athlete’s variation in RSI? What is a meaningful change in RSI for the team and for individual athletes?

# Some basic visualizations

VALD_ForceDecks_clean <- VALD_ForceDecks_clean |>
    filter(Date >= as.Date("2024-10-10"))
  
# Boxplot showing RSI distribution for each athlete
ggplot(VALD_ForceDecks_clean, aes(anon_id, RSI..m.per.s., color = anon_id)) +
  geom_boxplot() +
  labs(title = "RSI Variability by Athlete") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Line and scatter plot showing RSI trends over time for each athlete
# Includes a black line for the team average RSI to compare individual performance to
ggplot(VALD_ForceDecks_clean, aes(Date, Player.Average.RSI..m.per.s., color = anon_id)) +
  geom_point() + 
  geom_smooth(method = "lm", se=FALSE) +
  labs(title = "Player RSI Over Time with Team Average Line") +
  geom_line(aes(y = Team.Average.RSI..m.per.s.), 
            color = "black")

`geom_smooth()` using formula = 'y ~ x'

ggplot(data=VALD_ForceDecks_clean, aes(x=Date, y=Player.Average.RSI..m.per.s.)) +
  geom_point() +
  geom_smooth(method="lm")

#taking out all of the IDs in the data set
IDs <- unique(VALD_ForceDecks_clean$anon_id)

#for all unique IDs, calculate Kendall rank correlation, plot player RSI throughout season 
for(i in 1:16){
  RSI <- unique(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.)
  kendall_cor <- cor(1:length(RSI), RSI, method="kendall") #calculating Kendall correlation

  #plotting each player's RSI throughout the season with their correlation and rough model of the trend
  p <- ggplot(data = VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],], aes(x=Date, y=Player.Average.RSI..m.per.s.)) + 
    geom_point() +
    geom_smooth(method="lm") +
    labs(title = kendall_cor, subtitle = IDs[i])
  print(p)
}

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

#for all the players, calculate mean and variance of their RSI measurements and plot
for(i in 1:16){
  #for every player calculate mean and variance using time series assumptions and functions from the time series functions file
  mv <- sample_mean_and_variance(unique(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.))
  
  #plotting RSI measurements along with calculated mean and 1 standard deviation above and below mean for each player along with their measured RSIs throughout the season
  p <- ggplot(data = VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],], aes(x=Date, y=Player.Average.RSI..m.per.s.)) + 
    geom_point() +
    geom_hline(yintercept = mv[[1]]) +
    geom_hline(yintercept = mv[[1]] + sqrt(mv[[2]]), color="#CFB87C") +
    geom_hline(yintercept = mv[[1]] - sqrt(mv[[2]]), color="#CFB87C") +
    labs(title = "RSI Measurements Throughout Season with Mean and 1 Standard Deviation",subtitle = IDs[i])
  
  #printing out each plot
  print(p)
}

VALD_ForceDecks_clean$Player.Above.1.SD <- rep(NA,565)
VALD_ForceDecks_clean$Player.Below.1.SD <- rep(NA,565)
VALD_ForceDecks_clean$Player.Above.2.SD <- rep(NA,565)
VALD_ForceDecks_clean$Player.Below.2.SD <- rep(NA,565)

for(i in 1:16){
  mv <- sample_mean_and_variance(unique(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.))
  sdu1 <- mv[[1]] + sqrt(mv[[2]])
  sdl1 <- mv[[1]] - sqrt(mv[[2]])
  sdu2 <- mv[[1]] + (2*sqrt(mv[[2]]))
  sdl2 <- mv[[1]] - (2*sqrt(mv[[2]]))
  
#checking if RSI value is below 1 standard deviation for each player
  VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Below.1.SD <- ifelse((VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.<sdl1), 1, 0)
  
#checking if RSI value is above 1 standard deviation for each player
  VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Above.1.SD <- ifelse((VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.>sdu1), 1, 0)
  
#checking if RSI value is below 2 standard deviations for each player
  VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Below.2.SD <- ifelse((VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.<sdl2), 1, 0)
  
#checking if RSI value is above 2 standard deviations for each player
  VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Above.2.SD <- ifelse((VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Average.RSI..m.per.s.>sdu2), 1, 0)
}

#making a plot to show potentially significant RSI measurements throughout the season
ggplot(data=VALD_ForceDecks_clean, aes(x=Date, Player.Average.RSI..m.per.s.)) +
  geom_point(color = ifelse(VALD_ForceDecks_clean$Player.Above.2.SD == 1, "palegreen", ifelse(VALD_ForceDecks_clean$Player.Below.2.SD == 1, "red", ifelse(VALD_ForceDecks_clean$Player.Below.1.SD == 1, "red4", ifelse(VALD_ForceDecks_clean$Player.Above.1.SD == 1, "palegreen4", "grey80")))), alpha = 0.4) +
  labs(title = "Significant RSI Measurements Throughout the Season") 

This plot shows each player’s measured RSI throughout the season. Points that are colored gray are points where the player did not have a significant difference in their RSI from their mean RSI throughout the season. Dark green dots show RSI measurements that are outside of 1 standard deviation for each player and light green dots are RSI measurements that are outside 2 standard deviations for each player. Along the same lines, dark red dots shown points in which the player had an RSI outside of 1 standard deviation from their RSI distribution and light red dots show points in which the player’s measured RSI is outside of 2 standard deviations of their RSI distribution. One interesting thing about this plot is that there are no measurements that are below even 1 standard deviation for any player after February starts. If anything, there seems to be a higher concentration of RSI measurements that are at least above 1 standard deviation from the mean in the last 5 measurements of the season, after February starts than in the rest of the season.

#cycling through all of the players and plotting their RSI measurements throughout the season
for(i in 1:16){
  #plot the measured RSI for each player throughout the season
  #dark green means point is above 1 standard deviation from mean
#light green means point is above 2 standard deviations from mean
#dark red means point is below 1 standard deviation from mean
#light red means point is below 2 standard deviations from mean
p <- ggplot(data=VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],], aes(x=Date, Player.Average.RSI..m.per.s.)) +
  geom_point(color = ifelse(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Above.2.SD == 1, "palegreen", ifelse(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Below.2.SD == 1, "red", ifelse(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Below.1.SD == 1, "red4", ifelse(VALD_ForceDecks_clean[VALD_ForceDecks_clean$anon_id == IDs[i],]$Player.Above.1.SD == 1, "palegreen4", "grey70"))))) +
  ylim(0.7311404, 2.963925) +
  labs(title = "Measured RSI Values Throughout the Season and Their Significance")

print(p)
}